%% solveTridiag
%  Version 1
%  Author: Adeyinka Lesi
%  Date: 8/14/15

function sol = solveTridiag(A,B)
% solveTridiag performs the matrix division B/A*, where A* is the
% tridiagonalization of A
% A: Nx3 matrix
% B: Nx1 vector
% sol: Nx1 vector

% solve using tridiagonal matrix algorithm
aug = [A B];
len = length(B);
% forward sweep
for j = 1:(len-1)
    aug(j+1,[2 end]) = aug(j+1,[2 end]) - aug(j+1,1)/aug(j,2) ...
        * aug(j,[3 end]);
    aug(j+1,1) = 0;
end
% reverse sweep
for rj = 2:(len)
    aug(len+2-rj,end) = aug(len+2-rj,end)/aug(len+2-rj,2);
    aug(len+2-rj,2) = 1;
    aug(len+1-rj,end) = aug(len+1-rj,end) - aug(len+1-rj,3) ...
        * aug(len+2-rj,end);
    aug(len+1-rj,3) = 0;
end
aug(1,end) = aug(1,end)/aug(1,2);
aug(1,2) = 1;

sol = aug(:,end);